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SUMMARY 


CTJRSHL  is  a  high-precision, finite  element, 
for  tha  linear  static  stress  analysis  of  thin  elastic 
shells  of  quite  general  shape.  The  element  is  a  fully- 
conforming  displacement  element  of  curvilinear  tri¬ 
angular  form.  The  theoretical  basis  of  the  element  is 
presented  in  fairly  complete  detail,  and  some  aspects 
of  the  organization  of  the  computer  program  are 
discussed.  The  performance  of  CURSHL  is  examined 
with  regard  to  accuracy  of  numerical  integration, 
compulation  time,  and  representation  of  rigid-body 
modes.  t 
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1.0  Introduction 


CURSHL  is  the  acronym  both  of  a  high-precision  finite 
element  for  the  static  analysis  of  thin  shells  and  of  the 
associated  computer  program.  The  element  is  a  fully  conforming 
displacement  element  of  curvilinear  triangular  shape  v/hich 
utilizes  the  same  interpolation  functions  and  generalized 
displacements  as  previous  high-precision  elements  for  plate 
bending,  plane  stress,  shallow  shells  and  cylindrical  shells. 
Continuous  Kirchoff  constraints  are  used.  The  computer  program 
generates  the  stiffness  matrix  and  consistent  load  vector  for  a 
single  element.  Variable  thickness  and  elastic  modulus,  as  v/ell 
as  thermal  loads,  can  be  accommodated  by  the  program  but  at 
present  the  program  is  limited  to  linear  elastic  isotropic  material. 
The  program  provides  the  option  of  generating  a  stiffness  matrix 
based  on  either  Koiter-Sanders  shell  theory,  Donnell-Vlasov  shell 
theory,  or  shallow  shell  theory.  Except  for  a  requirement  of 
smoothness  the  shape  of  the  shell  is  virl  ually  unrestricted. 

A  general  description  of  CURSHL  was  first  given  at  CANCAM 
*71  [10]*.  An  amplified  description  and  an  outline  of  the 
considerations  which  guided  the  development  of  CURSHL  were  given 
in  Reference  [1],  which  also  included  several  examples  of  CURSHL' s 
performance.  Because  the  emphasis  in  these  references  was  on  the 
performance  and  accuracy  of  the  element,  many  details  of  the 
theoretical  and  computational  aspects  were  omitted.  The  purpose 
of  the  present  report  is  to  describe  in  fairly  complete  detail 
the  underlying  theory  and  to  outline  the  organization  of  the 
computer  program.  In  the  pages  which  follow  the  relevant  items 
from  shell  theory  are  presented,  formulas  for  the  interpolation  of 
displacement  within  an  element  are  derived,  and  the  method  of 
calculation  of  the  stiffness  matrix  is  described.  In  addition, 
certain  computational  details  of  the  program  are  discussed. 

Finally,  the  element  is  tested  in  a  small  number  of  example  problems 


* 


Numbers  in  square  brackets  denote  references  at  the  end  of  the 
text . 
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in  order  to  illustrate  its  accuracy. 

2.0  Rdsumd  of  Formulas  from  Shell  Theory 

The  required  formulas  from  shell  theory  are  assembled 
here  for  ease  of  reference.  General  tensor  notation  is  used  in 
order  to  avoid  unnecessary  restrictions  on  the  shape  of  ti.e  shell 
or  on  coordinate  systems.  Of  the  multitude  of  existing  first- 
order  shell  theories,  the  theory  of  Koiter-Sanders  can  claim., 
with  reasonable  justification,  to  be  the  best  [4].  It  is 
therefore  used  as  the  main  theoretical  basis  of  the  element.  It 
happens  that  the  equations  of  the  simplified  theories  of  Donnell- 
Vlasov  and  of  shallow  shells  are  obtainable  from  the  Koiter-Sanders 
equations  by  a  few  simple  modifications.  Since  the  latter  vneories 
are  widely  used,  they  are  made  optionally  available  in  the  computer 
program. 

The  treatise  of  Green  and  Zerna  [2]  Is  the  source  of  the 
formulas  for  the  geometric  parameters,  including  the  approximate 
formulas  which  may  be  used  when  the  shell  is  shallow.  The 
equations  of  Koiter-Sar.ders  shell  theory  are  taken  from  the  papers 
of  Koiter  [3],  and  of  Budiansky  and  Sanders [4'].  Tensor  notation 
following  the  conventions  of  Reference  [2]  is  used  throughout  this 
section. 

2 . 1  Geometric  Parameters 
( i )  General  Shell 

Let  0i,  82  be  curvilinear  coordinates  on  the  shell’s 
middle  surface.  Let  the  middle  surface  be  defined  by  the  position 
vector  r  =  r( 61,82)  whose  Cartesian  components  are  x,  ,  z.  The 
covariant  metric  tensor  is 

%  =  =  x-Xx-w  +  V’XV’V  *  Z’A%  (2-1) 

while  the  contravariant  metric  tensor  is  given  by 

a11  =a22/a,  a12  =  a21  =  -ai2/a,  a22  =  an/a  (2.2) 
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where 


a  -  a!i£22  ~  ai2a2i 


The  unit  normal  vector  is 

n  =  (r,i  x  rj2)/|r, i  x  r,2| 


and  we  note 


r,  i  x  r ,  2  |  =  /T 


The  covariant  and  mixed  curvature  tensors  are 


b>.p  =  *  V>Ap  +  V'Am  +  nz*’W 


A  Ap 
b  =  a  b 
y  PP 


(2,3) 


(2. 4) 


(2.3) 


(/■  .6) 


The  contravariant  base  vectors  are  given  by 

■+A  Ap-»- 

a  =  a  Kr- 


"P 


and  the  Christoffel  symbols  are 


rA  _  a^*r  =  a^pr  *r 
yv  Jyv  *p  syv 

=  aXP<x’px»Mv  +  ^’yv  +  z,pz,y\>^ 


(2.7) 


(2.8) 


The  formulas  for  covariant  differentiation  of  first  and  second 
rank  covariant  tensors  are 


AA|y  AA,y  rAyAp 
AAyjv  "  AAy,v  “  ^AvApy  ”  ^yvApA 


(2.9) 

(2.10) 


In  addition  to  the  foregoing  standard  formulas,  expressions 
for  the  covariant  derivatives  of  the  mixed  curvature  tensor  are 
needed.  To  obtain  these  we  first  differentiate  (2.6)  to  get 


-y  -y 


bMi,V  =  +  nVr’Av, 


(2.11) 


-  ij  - 


The  derivatives  of  the  normal  vector  are  given  by  the  Gauss- 
Weingarten  formula 


-*■  .  -c-y 

=  -V-' 


=  -a 


\>p 


so 


(2.12) 


Xp,v 


-y  -y 

n«  r 


1  Apv 


TO+  -*■  . 

2  PVV\>P 


-J* 

n*  r. 


Apv 


rp  D 
Ap  vp 


(2.13) 


Hence,  from  (2.10)  and  (2.13), 


bAp|v  n*ri,Apv  ^Ap^pv  ^Av^pp 


-  Tp  b 


pv  pA 


(2.14) 


It  is  interesting  to  note  that  b^|v  is  unaltered  by  any  permutation 
of  the  indices.  Since  the  covariant  derivatives  of  the  metric 
tensor  vanish,  the  required  formulas  are 


V|v 


_  „Xt. 
=  a  b 


Xp|v 


=  aAT(S.?, 


Apv 


-  it  b  -  rtr  b 


Ap  pv 


(2.15) 


-  Tp  b 


Av  pp  pv  pA 


) 


(ii)  Shallow  Shell 

When  the  shell  is  shallow  certain  approximations  are 
permissible,  as  discussed  in  Section  11.3  of  Reference  [2].  The 
metric  tensor  of  the  shell  can  be  approximated  by  the  metric 
tensor  of  the  base  plane,  and  hence  (2.1)  is  replaced  by 

a\i; =  x>xx>„ +  y’x5%  (2-16) 

The  contravariant  metric  tensor  is  still  given  by  (2.2)  but  now 
there  is  the  alternative  formula  for  a, 


a  =  (x,xy,2  -  x , 2 y , i ) 2 


(2.17) 
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Hence.,  within  the  shallow  shell  approximation, 

i:r,a  x  r,2l  =  /a  =  (x,iy,2  -  x,2y,!)  (2.13) 

where  the  correct  sign  of  the  expression  for  /a  was 
determined  from  the  particular  case  0j  =  x,02  =  y. 

The  unit  normal  can  still  be  calculated  from  (2.*0, 
and  in  view  of  (2.18'),  the  result  agrees  with  the  usual  shallow 
shell  approximation  for  the  normal.  That  is,  the  z-component 
of  n  is  1  while  the  other  two  components  are  equal  to  the 
negative  gradient  of  z.  If,  in  turn,  this  result  is  used  in 
(2.6)  to  calculate  the  curvature  tensor  the  expressions 
obtained  agree  with  the  usual  shallow  shell  approximation, 

bAp  *  zHam  (2-19) 


where  the  double  stroke  denotes  co variant  differentiation  j.n 
the  base  plane. 

The  Christoffel  symbols  for  the  base  plane  are  obtained 

jy  ->•  jj. 

by  replacing  r  in  (2.8)  by  R,  where  R  is  the  component  of  r 
lying  in  the  base  plane,  at  the  same  time  using  (2.16)  to 
calculate  the  metric  tensor.  Thus 


=  ^o’Hv  =  aXP(X>pX>MV 


+  y.py>IJU> 


(2.20) 


summary,  for  shallow  shells  formula  (2.1)  is 
replaced  by  (2.16)  and  formula  (2.8)  by  (2.20),  but  the  formulas 
for  a^,  n,  and  remain  valid. 


Membrane  and  Bending  Strains 


The  membrane  strain  tensor  is  related  to  the 
displacements  by 


:Ap  =  i(uA|p  +  um|A>  -  bApw 


(2- 21) 


=  V(uA,p  +  VP  '  rApUp  -  bApw 


o 


and  this  relation  applies  to  all  three  of  Koiter-Sanders , 
Donnell- Vlasov,  and  shallow  shell  theories. 

In  Koiter-Sanders  theory  the  rotations  of  the  normal 
and  the  surface  rotation  are  given  by 

^x  ‘  -  W’X  -  bXup 


(2.22) 


w,  =  r(u.  |  -  u  i,)  =  4(U-,  -  u  .) 

Ay  gv  A.|y  U I  A7  2V  A,y  y,A7 


(2.23) 


and  the  bending  strain  tensor  is 

k,  =  r(4»,  i  +  <*>  i , )  +  |(b?w  +  bpw,  ) 

Ay  2'Tx|y  y|A7  2  A  yp  y  Ap' 

=  -  w.,  +fp  w, 

’Ay  Ay  ’ p 


(2.2-4) 


-  ♦  bSup,x> +  +  b^x,p> 

+  -  b$|X  +  r?MbX+  r»p 


In  Donneil-Vlasov  theory  and  in  shallow  shell  theory  the 
rotations  of  the  normal  are  approximated  by 

=  -  w,.  (2.25) 


and  the  bending  strain  tensor  is  given  by  the  relatively 
simple  formulas 

KAy  -  2^A|y  +  ^y|A^  ”  “  wjAy 
=  -  w,Xy  +  rXMW>p 


(2.26) 


These  strain-displacement  relations  may  be  summarized  by  the 
matrix  relation 


(e)  =  [B]{d} 


(2.27) 
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where  {e}  is  a  column  'sector  of  strains, 

f*d 

{2}*  =  {s;  l,ei2jS223K:2i,iCi2,JC22-} 

and  {d}  is  a  column  vector  of  displacements  'md  their 
derivatives, 

{d}  =  {Ui,Ui,i,ar,2,il2.U2,l,U2,2,W,W,i,W,2vlW,Sl,W,J2,W,22} 

The  matrix  £8]  is  displayed  in  Table  1, 

2 , 3  Stress-Strain  delations  and  Strain  Energy 

The  stress  resultants  and  bending  moments  are  related 
to  the  membrane  and  bending  strains  ly 


am  =  chAmpt£ 

pt 

(2.30) 

>  =  DHXppT£ 

PT 

(2.31) 

where  C  is  the  stretching  rigidity,  D  the  flexural  rigidity, 
and 

HXppT  =  i«l-v)(aXVp  +  aXpapT)  +  2vaXpapT>  (2.32) 

For  homogeneous  shells  of  thickness  t, 

G  =  Et/(l-v?),  D  =  Et3/12(l-v2)  (2.33) 

For  laminated  or  sandwich  shells  other  expressions  for  C  and  D 
may  be  appropriate  and  their  use  is  not  ruled  out.  However 
(2.30),  (2.31),  and  (2.32)  apply  only  to  isotropic  shells  in 
which  there  is  no  coupling  in  the  stress-strain  relations 
between  bending  and  stretching. 

The  strain  energy  density  is 
dU/dfi  =  +  n>'"VAli) 

=  2(CHXppTe,  £  +  DHXppT 

z'  Am  pi 


(2.28) 


(2.29) 


k.  K  .. ) 
Am  p? 


(2.3*0 
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and  this  may  be  written  in  matrix  shorthand  as 

dU/dA  =  |{s}T[E]{e}  (2-.  35) 

The  matrix  [E]  is  displayed  in  Table  II. 

When  thermal  expansion  is  considered  the  s/  vess-s train 
relations  are  modified  to 

nXy  =  CHAypTepT  -  (Ea/(l-v))TRaAy  (2.36) 

mMi  =  j)HAyPTicp^  -  (Ea/(l-v))TMaAy  (2.37) 

where  a  is  the  coefficient  of  thermal  expansion  and  TR  and 
T^  are  the  resultant  and  the  moment  of  the  temperature  T  across 
the  thickness.  That  is 


|  T(c)dc 

-t/2 


(2.38) 


t/2 

Tm  =  j  ?T(5)d?  ’(2.39) 

-t/2 

where  5  is  the  coordinate  in  the  thickness  direction  of  the  shell. 
Formulas  (2.36)  and  (2.37)  apply  to  homogeneous  shells  and  would 
have  to  be  modified  for  laminated  and  sandwich  shells.-  The 
thermal  strains  do  not  affect  the  expression  for  strain  energy,  but 
instead  influence  the  expression  for  the  virtual  work  of  the 
external  loads. 

We  suppose  that  the  shell  is  acted  upon  by  a  distributed 
load  per  unit  area  p.  The  contravariant  components  of  p  are  pA 
in  the  tangent  plane  and  pn  normal  to  the  shell.  The  virtual 
work  per  unit  area  of  the  applied  loads,  including  the  contribution 
of  thermal  effects,  is 

dV/dA  =  p\  +  pnw  +  (Ecc/(l-v))(TRaAyeAM  +  T^e^)  <2.40} 
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or,  in  matrix  notation 

dV/aA  =  {P}T{d)  +  {Q}T{e>  (.2.41) 


where 

(P}T  =  0,0, 0,0,0}  (2.-42) 

{Q}T  =  (Ea/  (1-v) )  (T^a1 1 ,2T^a* 2  ,TRa22  jT^a1 1 ,2T^a1 2  ,T^a22  }  (2.43) 

2 . 4  Physical  Components 

The  tensor  components  of  stress  or  displacement  are  not 
the  same  as  the  physical  components  of  these  quancities.  It  is, 
of  course,  the  physical  components  which  ar^  ultimately  sought, 
although  tensor  component's  are  used  in  setting  up  the  stiffness 
matrix. 

The  physical  components  of  displacement  u,  v,  w,  resolved 
along,  the  covariant  base  vectors,  are  related  to  the  tensor 
components  by 

u  =  (alxui  +  alzu2)/ai i 

v  =  (a12Ui  +  a22u?. )  /a22  (2,4*0 


w  =  w 


The  contravuriant  tensor  components  of  load  are  related 
to  the  physical  components  pi,  p2,  pn,  by 


p1  =  Pi/v/aT7  ,  p2  =  p i/fsilz  ,  pn  «  5 


(2.45) 


The  physical  stress  resultants  and  the  physical 
bending  moments  are  related  to  their  tensor  counterparts 
by 
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Nn  =  n1  V(ai  1/a1  *) , 

N i 2  =  n12/a. 

N 2  2  =  n22/(a22/a22) 

(2.46) 

Mu  =  m11/(an/a11) , 

Mi  2  =  :m1?/a. 

M22  =  m22/(a22/a22 ) 

(2.4?) 

Tlis  physical  components  of  the  rotation  of  the  normal 
$i3$2,  are  related  to  the  tensor  components  by  the  formulas 


4>i  =  ,-a 1 1  <J>i  +  a12<f>2)/aii 
?2  -  (a21tj>i  +  a22<j>2)/a22 


(2.48) 


which  are  the  same  as  the  corresponding  formulas  for  the 
displacements.  According  to  Glockner  [12]  the  physical  component 
of  the  surface  rotation  ft  is  given  by 


» *  *“V/2 

=  (u2 , i  -  Ui , 2 )/2/a 
where  ea^  is  the  alternating  tensor. 

We  note  in  passing  that  the  rotation  vector  ft  is  given  by 


(2.49) 


3  =  ea(* 


ftn 


(2.50) 


The  conventions  for  the  positive  directions  of  the 
physical  components  of  displacement,  rotation,  stress  resultants, 
and  bending  moments  are  illustrated  in  Figures  2  and  3-  In  the 
figures,  the  surface  coordinates  have  been  denoted  as  a,  g', 
rather  than  8i,  02.,  The  positive  direction  of  the  normal  is  such 
that  the  normal  vector  forms  a  right-handed  triad  with  the 
tangent  vectors  along  the  a-curves  and  the  $-curves.  Figures 
2  and  3  are  intended  to  show  non-orthogonal  coordinates.  The 
directions  of  the  stress  resultants  are  parallel  to  the 
coordinate  curves.  The  vectors  which  represent  the  bending 
moments  are  orthogonal  to  the  coordinate  curves,  and,  consistent 
with  this,  the  components  of  bending  stress  are  parallel  to  the 
coordinate  curves. 
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3.0  Interpolation  Functions  for  Displacements 

At  this  -point  we  drop  the  tensor  notation  of  the  preceding 
sebtiori  and  introduce  a, 3  to  denote  the  general  curvilinear 
coordinates  on  the  shell  surface,  and  u,  v,  w  to  denote  the 
covariant  tensor  components  of  displacement.  The  coordinates 
a, 3  need  not  be  principal  nor  orthogonal. 

The  curvilinear  triangular  element  on  the  shell  surface 
can,  be  imagined  mapped  on  to  the  a-3  plane.  The  edges  of  the 
element  are  defined  by  specifying  that  they  are  straight  lines 
in  the  a-3  plane. 

The  displacements  in  CURSHL  are  interpolated  by  the  same 

functions  which  have  been  used  in  previous  high-precision 

elements  [5,6, 7, 8],  The  in-plane  displacements  u,v  are  taken 

as  complete  cubic  polynomials  in  a,  3,  while  w  is  taken  as  a 

restricted  quintic  polynomial  in  a, 3.  In  keeping  with  this  choice 

of  displacement  functions  the  generalized  displacements  are  the 

values  of  u,  3u/3a,  3u/33,  v,  3v/3a,  3v/33,  w,  3w/3a,  3w/33, 

32w/3a2,  32w/3a33,  32w/332  at  the  three  vertices  of  the  element, 

for  a  total  of  36  degrees  of  freedom.  Centroidal  displacements 

u  ,v  .  are  used  temporarily  as  degrees  of  freedom  during  the 
c  c 

development  of  the  stiffness  matrix  but  are  later  eliminated  by 
static  condensation. 

The  choice  of  interpolation  functions  and  generalized 
displacements  assures  that  the  element  is  fully  conforming 
provided  the  shell  surface  satisfies  certain  smoothness  conditions. 
For  Koiter-Sanders  theory  these  conditions  are  that  the  shell  must 
be  smooth  and  have  continuous  curvatures.  For  Donnell-Vlasov 
theory  and  for  shallow  shell  theory  it  is  sufficient  that  the  shell 
is  smooth.  These  conditions  are  discussed  in  more  detail  in 
References  [1]  and  [7]. 

The  order  of  the  discretization  error  of  a  conforming 
finite  element  depends  on  the  choice  of  interpolation  functions, 
and  can  be  predicted  by  the  Taylor’s  series  test  introduced  by 
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McLay  [11] ,  The  application  of  the  test  to  the  present  element 
is  essentially  the  same  as  in  tests  of  previous  high-precision 
elements  for  flat  plates  [5]  and  shallow  shells  [7].  The  result 
is  that  the  discretization  error  is  of  order  h6  where  h  is  a 
typical  linear  dimension  of  an  element.  The  Taylor" £  series  test 
also  shows  that  the  interpolation  functions  for  w  on  the  one  hand 
and  for  u,v  on  the  other  are  matched,  in  the  sense  that  they 
provide  equal  orders  of  accuracy  in  the  strains. 

Further  discussion  of  the  choice  of  interpolation 
functions  is  contained  in  Reference  [1]. 


3.1  Local  Coordinates 

Let  the  vertices  of  the  element  be  numbered  1,2,3  in 
counterclockwise  order  when  viewed  from  the  positive  side  of  the 
element.  Let  the  shell  coordinates  of  the  vertices  be 
(ai,3i),  (0(2,32)}  (013,03) i  respectively  as  shown  in  Figure  1. 

In  order  to  develop  formulas  relating  the  displacements 
within  an  element  to  the  generalized  displacements  it  is  helpful 
to  introduce  a  system  of  local  oblique  coordinates  £,n.  The 
S-axis  is  taken  to  coincide  with  the  1-2  edge  of  the  element  and 
the  n-axis  with  the  1-3  edge,  as  shown  in  Figure  1.  Further,  the 
local  coordinates  are  scaled  so  that  £  runs  from  0  to  1  along 
edge  1-2  while  n  nuns  from  0  to  1  along  edge  1-3.  The  relation 
between  a, 3  and  £,n  is  easily  established  as 


a  =  ai  +  (a2-ai  K  +  (a3-cti)n 
3  =  0i  +  (32-3iK  +  (3  3-3 1  )n 


(3.1) 


Derivatives  with  respect  to  local  and  global  coordinates 
are  related  by 
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3/H  =  f  i  i3/3a  +  f 2 i3/33 
3/3n  =  f  i  23/3(x  +  f223/33 

32/3C2  =  fn232/3a2  +  2fnf2132/3a33  +  f2i232/332  (3.2) 

32/3£3h  =  fnfi232/3a2  +  (fnf22  +  f  2  l  f  i  2 )  32/3a33  +.  f2if2232/332 
32/3n2  =  f  i2232/3ct2  +  2f  i2f2  23  2/3a33  +  f22232/332 

3/3a  =  dn3/3£  +  d2i3/'3n 

3/33  =  di 2  3/3£  +  d2  23/3n 

32/3ot2  =  dn  23  2/3?2  +  2dnd2132/3C3n  +  d2i  232/3n2  (3.3) 

32/3a33  =  dndi232/352  +  (dnd22  +  d2 id  1 2 ) 32/3£3n  +  d2id2232  /3n2 
32/332  =  d1  2  232/3?2  +  2di2d2232/3^3n  +  d2  2  23  2/3r)2 

where 


fll 

=  a2 

-  ot  1 

f  2  1 

-  3  2 

-  3 1 

di 1 

=  (3 

3  -  3 1 )/f 

d2  1 

=  _ 

(32  -  3i )/f 

f  = 

(a2 

-  ct  1 )  ( 3  3  - 

The  Jacobian  of  the  transformation 

3(a,3) 
3(5,n)'  = 


f  1 2  =  ct 3  -  aj 
f 2 2  -  33  -  3 1 

di2  =  -  (a3  -  oii)/f  (3**0 

d2  2  =  (.a 2  -  ai)/f 

1)  -  (a3  -  ai)(32  -  3 1 ) 
is  equal  to  f, 

f  (3.5) 
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3.2  Tangential  Displacements 

The  tangential  displacement  u  is  taken  as  a  complete 
cubic  polynomial  in  the  shell  coordinates  a, 3,  or  equivalently 
as  a  complete  cubic  polynomial  in  the  local  coordinates,  thus 

u  =  ai  +  a2£  +  a3n  +  +  as5n  +  a6n2  +  a?£3 

+  a8£2n  +  a9£n2  +  ai0n3 


(3.6) 


The  corresponding  generalized  displacements  are  the  values  of 

u  and  its  first  derivatives  with  respect  to  a, 3  at  the  three 

vertices  of  the  element.  In  order  to  match  the  number  of 

generalized  displacements  wi'.ti  the  number  of  coefficients  in 

(3.6)  the  displacement  u  at  the  centroid  is  introduced  as  an 

c 

additional  generalized  displacement.  Later  u  will  be  eliminated 

\y 

by  static  condensation. 

Denote  the  values  of  u,  9u/9g,  etc.  at  vertex  number  1  by 
Ui,  u^i,  etc.  Since  the  local  coordinates  of  vertices  1,2,3,  and 
of  the  centroid  are  (0,0),  (1,0),  (0,1),  (1/3, 1/3),  respectively, 
the  following  relations  are  obtained  from  (3.6) 

Ui  =  ai 


V  =  a2 

uni  =  a3 

u2  =  ai  +  a2  +  ai*  +  a? 

u^g  =  a2  +  2ae,  +  3a  7 

un2  =  a3  +  as  +  a8 

u3  =  ai  +  a3  +  a6  +  aio 


(3.7) 


u^3  =  a2  +  as  +  a9 
un3  =  a3  +  2a8  +  3aio 

u  =  ai  +  (a2+a3  )/3  +  (ai»+as+a6  )/9  +  (a7+a8+a9+ai o  )/27 

c 
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These  relations  may  be  explicitly  inverted  to  give 

{ Au>  =  [TU]{W^}  (3.8) 

where  {Au},  {W^}  are  the  column  vectors' 

T 

{  —  { 8. 1.^  3.2  j3-3^*»*3ioJ  (3*9.) 

{W^}T  =  {Ui>U^i}Uni,U2sU€2>Un2JU3.,U53JUn3,U0}  (3.10) 


and  the  matrix  [Tu]  is  displayed  in  Table  III. 

Let  (Wq)  be  the  column  vector  of  generalized  displacements 

{WU}T  =  {ui,uia,ui3,u2,u2a,u2g,u3,u;.a>u3g,uc}  (3.10a) 

where  Ui  ,  etc.  denotes  the  value  of  Su/8a  at  vertex  number  1, 
etc.  The  relation  between  (Wu)  and  {W^}  is  evident  from  equations 
(3.2.)  and  may  be  written  in  matrix  form  as 

{W&}  =  CRU]{WU}  (3.11) 

where  [Ru]  is  given  in  Table  III.  Hence  the  polynomial 
coefficients  of  (3.6)  are  related  to  the  generalized  displacements 
by 


{Au>  =  [TU][RU]{WU>  (3.12) 

Given  the  generalized  displacements,  the  polynomial  coefficients 
can  be  calculated  from  (3.12)  and,  in  turn,  the  values  of  u  and 
its  derivatives  can  be  found  from  (3.6).  A  similar  relation 
applies  to  the  displacement  v. 

An  advantage  of  using  local  coordinates  is  that  [Tu]  can 
be  found  explicitly  and  is  simple  in  form. 
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3.3  Normal  Displacement 

The  normal  displacement,  vr  is  taken  as  a  quintic 
polynomial., 

w  =  b!  4  b2£  +  b3n  b^2  +  ...  +  b 2 o Cn 4  +  b.21-n5  (3.13) 

As  in  previous  high-precision  elements,  constraints  are  imposed 
on  (3.13)  to  restrict  the  derivative  of  w  in  the  direction 
perpendicular  to  any  edge  to  a  cubic  variation  along  the  edge. 

The  notation  etc.  is-  used  to  denote  the  values 

of  w,  3w/3£.,  etc.  at  vertex  1,  and-  so  on.  The  following  -relations 
involving  twelve  of  the  polynomial  coefficients  are  obtained 
directly  from  (3.13) 


Wi  =  b 


"a  = 1)2 


wni  =  b> 


w55l  =  2b 4 


%1 =  bs 
wnni  =2bs 


(3.1*0 


w 2  =  b i  +  b2  +  bn  +  b?  +  b 1 1  +  bj 


w^g  =  b  2  +  2b 4  +  3b 7  +  4b ii  +  5b  is 

w^2  =  2b4  +  6b 7  +  12b  1 1  +  20bi6 

w3  =  bi  +  b3  +  bs  +  bio  +  bis  +  b2i 


w  ^  =  b 3  +  2b 6  +  3b 1 0  +  4b 1 5  +  5b 21 


w  ^  =  2b  6  +  6b  1 0  +  12b  is  +  20b  21 
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These  relations  can  be  explicitly  inverted  to  give 


b  i 
b2 

b  3 

b4 

bs 


=  w  l 
=  w 


SI 


=  w 


ni 

=  w£gi/2 


w 


Sni 


be 


(3.15) 


b7  =  -10Wi 
bio  =  -10wi 

b  1 1  =  15wi 
bis  =  15w i 
b  1 6  =  -  6wi 
b2J  =  -  6wi 


+  8w 


Si 


+  8w 


m 


-  ( 3/2)'^^  +  10W2  rr  +  W^/2 

-  <3/2)wnm  +  10Wi  -  ""ns  +  wnn3/2 

+  (3/2)w?w  -  15w,  +  7w52  -  wH2 

+  (3/2)wnm  -  15ws  +  7wn3  -  wnn3 

-  wK1/2  +  6w2  -  3w52  +  w552/2 

-  w  /2  +  6w3  -  3w  +  w  ,/2 

nni  n3  nn3 


To  obtain  formulas  for  the  other  coefficients,  the 
constraints  on  the  derivatives  of  w  must  be  considered.  Since 
terms  of  quartic  and  lower  degree  give  first  derivatives  which 
vary  at  most  cubically,  the  constraints  affect  only  the  six 
quintic  terms  in  (3.13). 

With  reference  to  Figure  1,  the  derivative  in  the 
direction  perpendicular  to  the  edge  1-2  in  the  a-$  plane  is 
given  by 


3w/3n  =  sin  0  3w/3a  -  cos  0  3w/33 


(3.16) 
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where 


sin  8  =  (82-3  x  )//s  3  3 >  cos  0  =  (a2-ai  ),.  Vs.3  3 
S33  =  (ci2-a  1 ) 2  +  (3  2-3  1 ) 2 

Transforming  (3.16)  to  derivatives  with  respect  to  £,n  ogives, 
after  some  manipulation 

3w/3n  =  {s 323w/3^-s 3 33w/3n}/(f/s3 3 ) 


(3.17) 


(3.18) 


where 


S32  =  (a2-ai ) (a3-ai )  +  (82-8 1 ) (8 3-8 1 )  (3.19) 

On  edge  1-2,  on  which  n=0,  the  contribution  to  the  term  in 
parentheses  in  (3.18)  from  the  quintic  terms  in  (3.13)  is  simply 

(5s 3 ?b x 6~s 3  3d  1 7 ) (3.20) 

and  hence  the  condition  for  cubic  variation  of  the  derivative 
perpendicular  to  the  edge  1-2  is 

bi7  =  (5s 32/S 3 3 )b 1 6  (3.21) 


The  above  equation  determines  bi?  since  bi6  is  already  known 
from  (3.15). 

Similarly,  the  constraint  along  the  edge  1-3  leads  to 
the  result 


b2o  =  (5s 3 2/s2 2 )b2 1  (3.22) 

where 

s22  =  (a.3-cxi)2  +  ( 8  3—8  x ) 2  (3.23) 

The  constraint  along  the  edge  2«3  leads  to  the  equation 

s 2 1  ( 5b 1 6-^b 1 7+3b 1 o-2b 1 a+b 2 0 )  =  s 3 j (b 1 7-au 1 $+3b 1 9-4b2o+5b2 1 )  (3.2*0 


I  >  #'**  /  If  f/l1  2, 


where 
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J-1"1 


S21  -  (a 3— ot  1 ) (a 3 — ct 2 )  +  (3 3-3 1 ) (3 3-32 ) 
s 3 1  =  (ct2-a  1 )  (a 3— ot2 )  +  (32-3 1 ) ( 3 3-82 ) 


(3.25) 


The  following  four  relations  are  obtained  directly  from 


(3.13) 


w 


n2 


w 


Sn2 


w 


53 


w 


Sn3 


=  b3  +  bs  +  b8  +  b  1 2  +  bi  7 


=  b  5  +  2b  8  +  3b  1 2‘  +  *lb  1 7 


=  b2  +  bs  +  b9  +  bi4  +  b2o 


=  bs  +•  2b  9  +  3b  1 1*  +  ^b2o 


(3.26), 


and  these  can  be  solved  to  give  bs,  bi2,  b9,  bi4  in  terms  of 
quantities  already  known. 


bs 

b  1 2  = 
b9  = 
b  1  4  = 


-  3V  -  2w5ni  +  3wn2  ‘  w5„2  +  b” 
2wni  +  w?m  "  2wn2  +  wen2  ‘  2017 

-  3wu  -  2v,s„i  +  3wu  ‘  wEn3  +  b;" 
2w5i  +  w5ni  "  2w?3  +  w5n3  "  2b2 ' 


(3.27) 


To  determine  the  remaining  coefficients  bi3,  bis,  bi9, 
there  are  the  relations  obtained  from  (3.13) 


w^r3  =  2(b4+b8+bj 3+bis ) 


wrin2  =  2(b6+b9+b 1 3+b 1 8  ) 


(3.28) 


and  the  constraint  equation  (3.2*0.  Solution  of  these  equations, 
making  use  of  (3.21)  and  (3.22)  gives  the  results 
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bi3  =  -  (2s2  i+3s  3 1  )w££g/2s  1 1  +  (3s2i+2s3i)wriri2/2sn 

+  (2S21  +  3S3l)('04+b8)/Sn  +  (  3S2  1+2S3  1  )  (bs+bg  )/S  1  1 

-  (Ai/sn)bi6  -  (A2/sii)bi7  (3.29) 

bie  =  w^^g/2  -  b &  -  b 9  -  bi3 
b 1 9  =  w.^/2  -  b 4  -  b8  -  b 1 3 

where 

sn  =  (a3-a2 ) 2  +  (3  3-^2 ) 2 

S21  =  (a3-ai )  (a3-ot2 )  +  (33-3 1 ) (3  3-32  ) 

S31  =  (a2-ai )  (a3-ot2 )  +  (32-3 1 )  (3 3-3z )  (3.30) 

A 1  =  -  5S21  +  5(4s.2l  +  S3l)S32/S33 

A2  =  -  5(S21  +  IIS3l)S32/S22  +  55  3  1 

In  obtaining  the  above  results,  the  fact  sn  =  s2i-s3i  was  used. 

The  connection  between  the  polynomial  coefficients  bi 
and  the  nodal  values  of  w  and  its  derivatives  can  be  expressed 
in  matrix  form 


where 


<v  ■ 


(3.31) 


{Aw}^  =  (b  1  ,b2 ,b 3 j . . . jb2 1 } 

{ww)T  =  !w‘-w5i’wni’w«i>w5ni>wnni’w:! 


(3.32) 
jWjj  <  n }  (3.33) 


In  the  computer  program  the  matrix  [T^J  is  set  up  in  the  following 
way.  First,  the  terms  in  twelve  of  the  rows  of  [T^ ]  are  implied 
by  the  relations  (3.15)  and  these  terms  are  inserted  into  the 
matrix.  In  view  of  the  relations  bi7  =  (5S32/S3 3)bi6  and 
b2o  =  (5S32/S22 )b2 1  the  rows  17  and  21  of  [T,r]  are  respectively 
equal  to  row  16  multiplied  by  (5s 3  2/s  93)  and  to  row  21 


multiplied  by  (5332/S22).  The  coefficient  b®  is  given  by  the 
first  of  formulas  (3-27),  that  is  by 


ba  -  -  3*ni  -  2W?I)1  +  3«„2  -  «en2  +  b,7 

Hence  to  set  up  row  8  of  [Tw3,  numbers  are  first  inserted  intvV 
this  row  corresponding  to  the  relation 


b  8 


+  3w 


H2 


W 


2 


and  then  row  17  is  added  to  row  8,  The  remaining  rows  9,12,13, 
14,18,19,  are  set  up  in  a  similar  way  using  relations  (3-27) 
and  (3-29). 

Let  (Ww)  be  the  column  vector  of  generalized  displacements 


lw„>  =  tw»wlo,wl8,www,o6 


,Wl30,v?2 


,w3: 


*} 


(3-34) 


where  wift  etc.  denotes  the  value  of  3w/3a  at  vertex  number  1. 
The  relation  between  (Ww)  and  {v/^}  is  evident  from  equations 
(3-2)  and  may  be  written  in  matrix  form  as 


(W'}  =  CRw](Ww} 


(3-35) 


where  [Rw]  is  displayed  in  Table  IV.  Hence  the  polynomial 
coefficients  of  v  5.13)  are  related  to  the  generalized  displacements 
by 


{Aw}  =  [Tw][Rw]{Ww)  (3.36) 

Although  [Tw]  has  not  been  written  down  explicitly,  its 
calculation  is  simple  and  requires  no  matrix  inversion. 
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4.0  Stiffness  Matrix 
4 . 1  Outline  of  Calculation 

The  expression  for  strain  energy  density  in  a  shell  was 
presented  in  section  2.2.  The  strain  energy  of  the  finite 
element  therefore  is 

Ue  =  £/7{e}T[ E]{e}/a  dad3  (4.1) 

where  /a  dadB  is  the  differential  of  area  on  the  shell  surface. 

The  terms  of  the  stiffness  matrix  [Kl  can  be  obtained  from  (4.1) 
by  noting  that  the  i-jth  term  of  [K]  is  given  by 

k, ,  =  32Ue/3X,.3X.  (4.2) 

■LJ  J 

where  denotes  the  i-th  generalized  displacement.  With  the 
notation 


.{ei>  =  3{e}/3Xi  (4.3) 

differentiation  of  (4,1)  gives,  in  view  of  the  symmetry  of  [E], 

kij  =  dad3  (4.4) 

Since  the  strains  depend  linearly  on  the  displacements  u,v,w, 
which  in  turn  are  linearly  proportional  to  the  generalized 
displacements  X^ ,  the  quantity  { e^ }  may  be  interpreted  as  the 
strains  which  occur  when  the  i-th  generalized  displacement 
has  unit  value  and  all  other  generalized  displacements  are  zero. 

The  evaluation  of  {e-^}  will  be  discussed  later. 

When  the  integration  is  transformed  to  local  coordinates 
formula  (4.4)  becomes 

=  //{ei}T[E]{e^  }f/a  dad3 


(4.5) 
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where  f  is  the  Jacobian  of  the  transformation.  The  strains  which 
appear  in  the  integrand  of  (4.5)  depend  in  a  complex  way  on  the 
shell's  curvature,  its  metric  tensor,  and  the  Christoffel  symbols,. 
The  latter  quantities  can  be  complicated  functions  of  the  shell 
coordinates.  Rather. than  attempt  to  approximate  them  by 
polynomials,  thus  making  closed-form  evaluation  of  (-4.5)  possible, 
it  seems  more  convenient  to  evaluate  (4,5)  by  numerical  integration. 

The  use  of  numerical  integration  introduces  a  certain  error. 
In  order  not  to  degrade  the  accuracy  of  the  high-precision  element., 
the  error  in  the  numerical  integration  formula  should  not  be 
worse  than  h6,  which  is  the  order  of  the  discretization  error. 

In  fact,  a  formula  with  an  error  of  order  h8  is  used,  in  an 
attempt  to  make  the  error  of  numerical  integration  negligible. 
Details  of  the  formula,  which  is  a  13-point  quadrature  formula  of 
Gaussian  type,  are  given  in  Table  V  and  Reference  [93.  The 
formula  is  fully  symmetric  with  respect  to  the  three  vertices. 
Application  of  the  formula  is  facilitated  by  the  fact  that  the 
local  coordinates  £,n  are  two  of  the  thi'ee  area  coordinates  of 
points  within  the  triangle  while  the  third  area  coordinate  is 
i-5-n. 

When  numerical  integration  is  applied  to  (4.5)  the 
formula  becomes 


=  1  1%  <<eilTCE]{eJ>r*/S)n  C.6) 

where  the  subscript  n  on  the  bracketted  quantity  indicates 
that  the  quantity  is  to  be  evaluated  at  the  n-th  pivotal  point 
of  the  numerical  integration.  It  is  this  formula  which  is  the 
basis  of  the  calculation  of  the  stiffness  matrix.  The  sequence 
of  calculation  is  as  follows.  At  the  first  pivotal  point  of  the 
numerical  integration  the  geometric  parameters  of  the  shell  are 
computed  and  the  column  vector  (e, }  is  evaluated  for  each  i  and 
stored.  Then  the  expression  |{ei>  [E]{ej}f/a  is  computed  for  all 
i  and  j ,  is  multiplied  by  the  weighting  factor  cn  of  the  numerical 
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integration 3  and  the  results  are  stored  in  the  array  which  the 
completed  stiffness  matrix  will  eventually  occupy.  The  procedure 
is  then  repeated  at  the  eond  pivotal  point  of  the  numerical 
integration  and  the  results  are  added  term-by-term  to  the 
quantities  which  already  occupy  the  stiffness  matrix  array.  The 
procedure  is  then  repeated  at  the  third  pivotal  point  ;  nd  continues 
until  all  pivotal  points  nave  been  covered.  In  this  way  the  uncon¬ 
densed  stiffness  matrix  is  built  up.  The  final  step  is  to 
statically  condense  the  stiffness  matrix  so  as  to  eliminate  the 
centroidal  displacements. 

The  foregoing  sequence  of  calculations  seems  to  be  as 
economical  as  any,  both  in  regard  to  storage  space  and  number  of 
operations . 

4 .  2  Evaluation  ofle^} 

It  was  shown  in  section  2.2  that  {e}  is  related  to  the 
vector  pf  displacements  and  their  derivatives  { d >  by 

{e}  =  [SHd}  (4.7) 


where,  in  the  notation  of  the  present  section 


f 

\ 


(u,u 


a 


a 


w  ,wD,w  ,w  D,W,.J 
*  3  aa*  a3  33 


a 


(4.8) 


the  subscripts  denoting  derivatives  with  respect  to  a  and  3. 
We  define  the  corresponding  vector  {d! }  involving  derivatives 
with  respect  to  £,  and  a. 


{d-}T  =  {uJu5JunJv,vE,vn,w,w?,Wn,w55>w5n,wiln} 


(1.9) 


From  equations  (3.3)  the  relation  between  {d}  and  (d* }  is 


{d}  =  [Rd]{d'} 


(4.10) 


where  [Rd]  is  given  in  Table  VI.  Hence 
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{e}  =  [B3[Rd]{d'}  (4.11) 

and  therefore 

{e1>  =  [B][Rd]{d|>  (4.12) 

where  {d^}  is  defined  as 

(d'}  =  3{d}/3X.  (4.13) 

Since  the  displacements  are  linearly  proportional  to  the 
generalized  displacements,  (dp  can  be  interpreted  as  the  value 
of  {d1}  when  the  i-th  generalized  displacement  has  unit  value  and 
the  other  generalized  displacements  are  zero. 

The  procedure  for  finding  (ei }  then  is  the  following. 

The  i-th  generalized  displacement  is  given  unit  value  and  all 
other  generalized  displacements  are  set  equal  to  zero.  The 
corresponding  polynomial  coefficients  are  calculated,  from  (3  12) 
and  (3..36).  The  values  of  displacement  and  their  derivatives, 
other  words  the  value  of  {d^} ,  are  then  calculated  from  the 
polynomial  interpolation  functions.  Finally  {e.}  is  obtained 
from  (4.12). 

4 . 3  Input  Data  for  Matrices  [B]  and  [E] 

The  matrices  [B]  and  [E]  depend  in  a  very  complex  way  on 
the  geometric  parameters  of  the  shell.  In  general  these 
parameters  vary  over  the  shell  surface,  and  hence  [B]  and  [E]  must 
be  evaluated  at  each  pivotal  point  of  the  numerical  integration. 

As  assumed  in  section  2.1,  the  shape  of  the  shell  is 
specified  by  giving  the  position  vector  ?  as  a  function  of  the 
shell  coordinates 


r  =  r  (a, 6) 
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or,  in  component  form,  by  specifying 

x  =  x(a,3) ,  y  =  y (a,3) ,  z  =  z(a,3) 

where  x,y,z,  are  the  Cartesian  components  of  r.  By  using  the 
formulas  of  section  2.1  all  terms  of  [B]  can  be  computed  if  the 
values  of  x,y,z  and  their  derivatives  are  known.  Derivatives  up 
to  third  order  are  required  in  Koiter-Sanders  theory  while 
second  order  derivatives  suffice  for  Donne 11- Vlasov  and  shallow 
shell  theory.  This  data  is  fed  into  the  computer  program  from  a 
user-supplied  subroutine  which  must  return  the  values  of  x,y,z, 
and  their  derivatives  at  an  arbitrarily  gj.ven  point  cc,3.  If,  as 
is  generally  the  case,  the  shell  surface  is  of  simple  form  then 
the  exact  equations  of  the  surface  can  be  used  in  setting  up  the 
subroutine,,  On  the  other  hand  the  representation  of  the  shell  by 
a  fitted  polynomial  surface  is  not  precluded. 

The  computer  program  incorporates  switches  to  evaluate 
[B]  according  to  either  of  the  three  shell  theories. 

The  matrix  [E]  depends  not  only  on  the  geometric 
parameters  but  also  on  the  elastic  parameters  C,D,  and  v.  Since 
numerical  integration  is  used,  variable  thickness  and  elastic 
moduli  can  be  handled  easily.  It  is  only  necessary  to  use  the 
appropriate  values  of  C,D,  and  v  at  each  pivotal  point  of  the 
numerical  integration.  This  data  is  fed  into  the  computer  program 
from  a  second  user-supplied  subroutine  which  must  return  the 
values  of  C,D,  and  v  at  any  arbitrarily  given  point  a, 3.  This 
method  of  handling  the  input  data  accommodates  variable  elastic 
properties  but  still  is  fairly  simple  in  the  case  of  constant 
properties . 
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4.4  Efficient  Calculation  of  Strain  Energy  Density 


In  the  computer  program  the  calculation  of  the  expression 
-^cn ( (ei }  [E]{ej )f»/a)n  in  (4.6)  occurs  within  a  triple  do-loop, 

.the  loops  spanning  the  rows  of  the  stiffness  matrix,  the  columns 
of  the  stiffness  matrix,  and  the  pivotal  points  of  the  numerical 
integration.  Consequently  each  operation  performed  in  evaluating 
the  expression  is  repeated  some  9600  times.  It  is  therefore  good 
-to  reduce  the  number  of  operations  to  a  minimum. 

If  the  expression  is  multiplied  out  it  is  found  that  it 
can  be  reduced  to 


icn{ei}T[E](ej}f/a 

=  ifcnC(eie'5v'a  +  ( (l-v)//a)  (2e\ 2**ei  ^  2“e22e'!  x  ^  (4.14) 

+  |fCnd(xiKJv^L  +  ((1-V)/Za)  (2K*2K*j2-K*il42  K22Kll)) 


where 


e1  =  a11^  +  2alze]z  +  a22e22 
k1  =  au<},  +  2a1  zk^2  +  a22^ 


(4.15) 


and  where  the  superscript  i  denotes  the  components  of  {e^}.  That 
is 

{e^}  =  (s j j , e j 2 , e 2 2 j j , k j 2 ,k 2 2 )  (4.16) 

If  we  define 

eh  a  A (l-v)cnCf/2/a)el 1 
£"22  =  /"( (1-v ) c^Cf/2/a) e 2  2 

i  i  (2j>17) 

i u  =  A (l-v)cnCf//a)ei2 

e1  =  /(cnCf/a/2)ei 
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Kll 


<2  2 


<\z 


K1 


then  expression  (4.14) 


fc^e/tEHe^f/a 


=  /■((l-v)cnDf/2/a)K}1 
=  A(l-v)cnDf/2/a)42 

=  /((l  -v)cnDf//a)icl2 
=  /(cnDf  /a/2)K* 

is  reduced  to 

=  +  E’ly.E'lz  “  clielz  **  eIzE^  1 

+  K1^  +  ic{2K^2  -  Kiik^2  -  ic22ic^r 


(4.18) 


(4.19) 


The  right-hand  side  of  (4.19)  requires  only  eight  multiplications 
and  is  the  formula  actually  used  in  the  program.  The 
multiplications  involved  in  forming  the  barred  quantities  in 
(4.17)  and  (4.18)  can  be  taken  outside  two  of  the  do-loops  and 
do  not  add  significantly  to  the  computation  time. 

Computing  time  is  also  saved  by  exploiting  the  fact  that 
the  bending  strains  do  not  depend  on  u  or  v  in  Donnell-Vlasov 
and  shallow  shell  theories.  Hence  in  these  cases  the  bending 
strains  can  be  omitted  from  (4.19)  whenever  the  i-th  or  j-th 
generalized  displacement  relates  to  u  or  v. 


5.0  Load  Vector 

The  calculation  of  the  load  vector  is  quite  similar  to 
the  calculation  of  the  stiffness  matrix.  The  formula  for  the 
virtual  work  per  unit  area  of  the  applied  loads,  including  thermal 
effects,  has  been  given  in  section  2.3.  The  total  virtual  work  of 
the  loads  acting  on  a  finite  element  then  is 


Ve  =  // ( {P}T{d>  +  {Q}T{e})/adadg 

=  //({P}TCRd]{d’ }  +  {Q}^{e} ) /adadg 


-  29  - 


Hence  the  i-th  term  of  the  load  vector  {L}  is  given  by 


i±  -  3Ve/ax± 

=  //({?}TL’RdKd|}  +  {Q}T{ei>)/adadB 


(5.2) 


After  transformation  to  local  coordinates  and  evaluation  by 
numerical  integration,  (5.2)  becomes 


*  1  l  cn((lP>TCRd){di>  +  {Q}T{ei})f/i)n  (5.3) 


To  compute  the  above  expression,  values  of  the  applied 
loads  and  temperatures  must  be  .known  at  each  pivotal  point  of 
the  numerical  integration.  This  data  is  fed  into  the  program  from 
a  user-supplied  subroutine  which  must  return  the  values  of  the 
three  physical  components  of  load,  together  with  the  temperature 
resultant  and  moment,  at  an  arbitrarily  given  point  a, 3.  For 
convenience,  this  subroutine  and  the  subroutine  which  furnishes 
the  values  of  the  elastic  constants  are  combined  into  one. 

The  technique  which  was  used  to  make  the  calculation  of 
the  stillness  matrix  more  efficient  is  also  used  to  simplify  (5.3). 

6.0  Performance  of  CURSHL 

In  this  section  some  results  which  illustrate  the 
performance  and  accuracy  of  CURSHL  are  presented.  It  goes  without 
saying  that  numerous  unreported  tests  were  also  run  with  CURSHL  to 
check  that  the  program  functions  correctly,  and  to  discover  and 
eliminate  bugs. 

6.1  Accuracy  of  Numerical  Integration 

In  order  to  obtain  some  appreciation  of  the  errors  due  to 
numerical  integration  the  problem  of  a  pressurized  spherical  cap 
was  solved  and  compared  with  a  previous  solution  based  on  the 
high-precision  shallow  shell  element  of  Reference  [7].  Both 
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elements  use  the  same  interpolation  functions  but  the  latter 
element  uses  exact  closed-form  integration,  so  the  discrepancies 
betvjeen  the  two  solutions  are  due  entirely  to  errors  of  numerical 
integration  in  CURSHL.  These  discrepancies  are  tabulated  in 
column  B  of  Table  VIJ.;  For  comparison,  the  error  in  the  finite 
element  solution  of  Reference  [7]  is  tabulated  in  columns  A  of 
the  table.  Thus  columns  A  represent  the  discretization  error 
of  the  high-precision^  elements  and  columns  B  the  additional  error 
in  CURSHL  due  to  numerical  integration.  In  general,  the  latter 
error  is  small  compared  to  the  former,  as  expected. 

The  geometry  of  the  problem  is  illustrated  in  Figure 

6 . 2  Computation  Time 

The  time  required  by  CURSHL  to  compute  one  stiffness 
matrix  and  load  vector  lies  between  1.1  and  1.7  seconds  depending 
on  the  nature  of  the  shell  and  the  choice  of  shell  theory.  More 
details  are  given  in  Table  VIII.  The  times  quoted  were  obtained 
on  an  IBM  360/67  machine  under  TSS  (Time  Sharing  System).  The 
programs  had  been  compiled  by  an  optimizing  compiler  using  H 
level  Fortran. 

6. 3  Rigid  Body  Modes 

Ideally  the  stiffness  matrix  should  have  six  zero 
eigenvalues,  corresponding  to  the  six  rigid-body  modes  of  the 
element.  The  smallness  of  the  first  six  eigenvalues  is  often 
used  as  a  measure  of  the  quality  of  the  stiffness  matrix.  The 
eigenvalues  of  CURSHL  for  a  number  of  cases  are  presented  in 
Table  IX. 

The  first  case  is  a  flat  plate,  and  the  effect  of  changes 
in  the  thickness  is  examined.  It  is  seen  that  the  first  six 
eigenvalues  are  satisfactorily  small.  Next  are  considered  three 
cases  where  the  strain-free  modes  can  be  represented  by  polynomial 
interpolation  functions  and  hence  should  be  exactly  represented 
by  CURSHL.  The  cases  include  a  flat  plate,  a  shallow  spherical 
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shell,  and  a  deep  cylindrical  shell  using  Donnell-Vlasov  theory. 
The  first  :six  eigenvalues  are  againv  satisfactorily  small  in  these 
cases.  Finally  cylindrical  and  spherical  elements  are  considered 
using  Koiter-Sanders  deep  shell  theory,  and  the  effect  of  changes 
in  curvature  is  examined.  In  these  cases  the  rigid  body  modes 
can  only  be  approximated  by  polynomial  interpolation  functions. 
Here  the  first  six  eigenvalues  are  reasonably  small  if  the 
curvature  of  the  element  is  shallow  but  become  larger  as  the 
curvature  increases. 

The  eigenvalues  quoted  here  differ  drastically  from  the 
eigenvalues  quoted  in  Reference  [1].  The  reason  is  that 
Reference  [1]  gives  the  eigenvalues  of  KX  =  AX,  K  being  the 
stiffness  matrix,  while  Table  IX  gives  the  eigenvalues  of 
KX  =  AMX,  where  M  is  the  consistent  mass  matrix.  The  computation 
of  M  follows  the  same  general  lines  as  the  computation  of  K  except 
that  a  36-point  numerical  integration  formula  was  used  in  this 
instance..  The  equation  KX  =  AMX  is  preferred  to  KX  =  AX  for  two 
reasons.  With  the  former  equation  the  ratios  of  successive 
eigenvalues  are  independent  of  scale,  and  the  eigenvalues  have  a 
physical  interpretation  in  terms  of  natural  frequencies. 

For  all  calculations  the  stretching  rigidity  Et/(l-v2) 
and  the  surface  mass  density  pt  were  taken  as  unity,  while  the 
form  of  the  element  was  a  right-angled  isosceles  triangle  with 
unit  sides. 

It  may  be  noted  from  Table  IX  that  the  first  eigenvalue, 
although  always  small,  is  sometimes  negative.  Therefore  the 
stiffness  matrix  does  not  quite  have  the  expected  property  of 
positive  semi-definiteness.  The  negative  eigenvalue  is  felt  to 
be  due  to  errors  of  numerical  integration  and  possibly  to 
round-off. 
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Table  I 

Strain-Displacement  Matrix  [B] 


In  Donne 11- Vlas ov  and  shallow  shell  theory 


Table  II 

Elasticity  Matrix  [E] 
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O  DEj 


[E,]  =  (a11)2 

2a1 1a12 


2as ia12  (l-v)(a12)2+va11^ 

2(l-v)a1 la22+2(3+v) (a12/2  2a22a12 
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Table  III 

Transformation  Matrices  for  Tangential  Displacements 
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Table  IV 

Transformation  Matrices  for  Normal  Displacement 

R2  0  0 

0  R2  0 

0  0  r2 


[R2]  =  ‘  1 

0 
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’  0 
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Numerical  Integration  Formula  Used  in  CURSHI 
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Table  VI 

Transformation  Matrix  for  Generalized  Displacements 
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Table  VIII 

Computing  Times  for  CURSHL 


Shape 

Theory 

VMR  time,  sec. 

flat  plate 

shallow 

1.1 

flat  plate 

K-S 

1.2 

circular  cylinder 

D-V  , 

1.1 

circular  cylinder 

K-S 

1.4 

sphere 

D-V 

1.3 

sphere 

K-S 

1.6 

elliptic  cylinder 

K-S 

1.4 

cone 

K-S 

1.5 

Table  IX 


Eigenvalues  of  KX  =  XMX,  where  K  is  the  CURSHL  stiffness  matrix 
(a)  Plat  plate  element;  effect  of  thickness 


EIGENVALUE 

■SB 

t  =  0.02 

Xi 

-o.7ixlo~14 

-0.70xl0“14 

-0. 70x10"" 1 4 

Xs 

0.70xl0“l" 

0.71xl0"14 

0.71xl0"14 

X7 

0.13X101 

0.1 3x10“ 1 

0.13xl0""3 

X36 

0.15*104 

0.32xl03 

0.32xl03 

X6/X7 

0.56x10"14 

0.36xl0"12 

0.56xl0“10 

(b)  Cases  where  exact  zero  eigenvalues  are  expected,  t=0.02;  R=1 


EIGENVALUE 

PLAT  PLATE 

SHALLOW 

SPHERE 

CYLINDER, 

D-V  THEORY 

Xi 

-0.70X10-14 

-O.llxlO""13 

-0.76X10'14 

Xs 

0.71X10"14 

0.80xl0“14 

0.77X10"14 

X7 

0.13X10""1 

0.l4xio_1 

0.14X10""1 

X3 1» 

0.32xl03 

0.4lxl03 

0.32x10s 

X6/X7 

0.56xl0"12 

0.59X10"12 

0.55X10""12 

. . . (cont 'd) 
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FIG.  4  DIAGRAM  OF  SPHERICAL  CAP.  UNIFORM  APPLIED 
PRESSURE.  ALL  EDGES  SIMPLY  SUPPORTED 
Rt/L2  =  0.02 


